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I discuss a family of statistical-mechanics models in which (some classes of) elements of a finite group Q 
occupy the (directed) edges of a lattice; the product around any plaquette is constrained to be the group identity 
e. Such a model may possess topological order, i.e. its equilibrium ensemble has distinct, symmetry-related 
thermodynamic components that cannot be distinguished by any local order parameter. In particular, if Q is a 
non-abelian group, the topological order may be non-abelian. Criteria are given for the viability of particular 
models, in particular for Monte Carlo updates. 

PACS numbers: 75.10.Hk,05.50.+q, 2.20.Hj 



I. INTRODUCTION 



"Topological order" llj,|2|] in a system means it has an emer- 
gent ground state degeneracy (in the thermodynamic limit), 
but (in contrast to symmetry -breaking), no local order param- 
eter operator can distinguish the states. Topological order has 
attracted great interest over the last 20 years, since (i) it can- 
not (by definition) be captured by the Landau order-parameter 
paradigm and is hence exotic from the viewpoint of traditional 
solid-state theory; J2l; (ii) it is associated with "fractional - 
ized" excitations; (iii) it is proposed to implement qubits by 
the ground-state degeneracy, the coherence of which is robust 
against environmental perturbations |H101; (iv) the formula- 
tion in terms of ground-state degeneracy makes it attractive 
for numerical exploration by exact diagonalization @1. 

The best-known examples of topological order are 
quantum-mechanical: the quantum Hall fluids) and lattice 
models based on Z2, the simplest group, such as Kitaev's toric 
code 0] ■ Indeed, Wen JH 0] once emphasized quantum me- 
chanics as a defining property of topological order. But we can 
separate these notions: topological order (as defined above) is 
meaningful in a purely classical model (as developed in this 
paper) or in a quantum-mechanical model at T > [9] (so that 
e.g. its renormalization-group fixed points represent classical 
behaviors). Indeed, I would suggest that the subject of topo- 
logical order skipped over more elementary examples, owing 
to historical accident. Compare with the history of sponta- 
neous symmetry breaking: theorists understood classical crit- 
icality long before quantum criticality II 1011 . and we approach 
quantum critical properties in light of their similarities or dif- 
ferences from the classical case. 

Analogously, it is hoped that, in the case of topological or- 
der, classical models will (at the least) be a pedagogical aid, 
and that behaviors evidenced in classical models may provide 
a framework for conjectures about the quantum models. Dis- 
entangling classical notions and inherently quantum mechan- 
ical ones might lead to clearer (or at least different) thinking. 
Also, the framework in this paper naturally draws us to face 
hitherto unfamiliar groups - e.g. the group A$ (see Sec. XX) 
- and it might inspire the construction of quantum-mechanical 
models involving these groups. 

The explicit notion of "classical topological order" was in- 
troduced and highlighted in ]9J], in particular the points that 
(i) it is characterized by ergodicity breaking, (ii) can be im- 



plemented by hard constraints, and (iii) must have a discrete 
classical dynamics - all of which applies to the models in this 
paper. However, much of Ref. U9[] was framed in terms of 
the relation of the classical model to a quantum model, e.g. 
by taking the quantum model to a temperature at which quan- 
tum coherence is no longer important, [40] or by removing 
some Hamiltonian terms. In the present work, the model is 
formulated from the start as an ensemble of classical statisti- 
cal mechanics, without concern for the existence of a quantum 
counterpart. That will permit consideration of a richer set of 
models (i.e. discrete non-Abelian groups Q), for which we 
might not know how to concoct a good quantized version. 

I will consider models based on either abelian or non- 
abelian discrete groups. It should be noted that abelianness 
in this paper has a different significance than in the quantum 
context. In the latter case, the group in question is the Berry 
phase (or its generalization to a unitary matrix) induced by 
evolution of the wavefunction from one state to an equivalent 
one. A specific case is the statistics of quasiparticles whose 
world lines braid around each other. Quantum-mechanical 
non-abelianness may be realized in models based on abelian 
groups such as Z%. Most of the fractional quantum Hall flu- 
ids are abelian, but non-abelian statistics is more suitable for 
quantum computation Proposed realizations of non- 

abelian topological order in this sense are formulated in lattice 
models as sums over loop coverings 

BUS. 

The primary focus here is not on ideas that point to analytic 
solutions or to connections with the existing literature of topo- 
logical order. Instead, this is meant as a generic blueprint for 
numerical studies. For example, the classification of models 
in Sec. [Hi] (as summarized in the tables) is motivated by the 
need to select a good one for simulations, and the quantities 
defined in Sec. [V] are all measurable in Monte Carlo simla- 
tions. However, the actual simulation results will be left to 
subsequent papers 111 111 . 



A. Height models 

I shall realize topological order by generalizing the concept 
of "height models". Their defining property IIL?T - [20ll is the 
existence of a mapping directly from any allowed spin con- 
figuration {/Ji} to a configuration of heights h(r), wherein 
h(r) — h(r'), for adjacent sites r and r', is a function of the 
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spin variables in the neighborhood (normally, either two spins 
on the sites r and r', or else on one spin on site i at the mid- 
point of the r-r' bond: spins and heights commonly live on 
different lattices). 

The "spins" in the model could be any discrete degree of 
freedom - e.g. dimer coverings. For {h(r)} to be well de- 
fined, it is necessary (and sufficient) that the sum of height 
differences is zero round any allowed plaquette configuration. 
Thus, a (local) spin-constraint is assumed that excludes (at 
least) the configurations without well-defined heights, yet still 
allows a nonzero entropy of states So in the thermodynamic 
limit. In the simplest cases, h(r) is integer-valued. 

Such models may have "rough" phases, in which the 
coarse-grained h(r) behaves as a Gaussian free field, i.e. the 
effective free energy of long-wavelength gradients is 




(1.1) 



Via the apparatus of the Coulomb gas formalism ll2lll22Tl . this 
implies the spin variables have power-law correlations (criti- 
cal state); the topological defects may have unbinding tran- 
sitions like the Kosterlitz-Thouless transition. Indeed, in two 
dimensions most critical states can be addressed as "height 
models" 112 If! and this formalism provides an alternative route 
to computing exact critical exponents 11711 . besides conformal 
field theory. 

Note that on a topologically non-trivial space (such as the 
torus periodic boundary conditions), there are nontrivial loops 
I, such that the net height difference (or "winding number") 
we added around such a loop is nonzero. It is easy to see this 
is a topological invariant, in that we. is unchanged if the loop is 
shifted and deformed (so long as it stays topologically equiva- 
lent to the original one.) Thus, if ... are the fundamental 
loops, the configuration space divides up into sectors labeled 
by (wg 1 , W£ 2 ,...). Here, and also for the discrete-group height 
models introduced in the paper, a sector is each set of config- 
urations which can be connected to each other by a succession 
of local spin changes (i.e. "updates" in the terminology used 
later). 

Point defects may also be admitted and loops around them 
may also have a nontrivial wi (in which case they are topo- 
logical defects). Clearly, the winding number wg in a height 
model is analogous to t € T in a topologically ordered model; 
we could almost say this is a special case of topological order 
in which T is Z, here meaning the (infinite discrete) group of 
integers under addition. 

In place of a Landau order parameter, the (near) degenerate 
states in a topologically ordered system may instead be distin- 
guished by a global loop operator, acting around a topologi- 
cally nontrivial loop I. Just as a Landau order parameter forms 
a group representation of a broken symmetry, and labels the 
symmetry-broken states, in topological order the global loop 
operator ought to form a faithful representation of the "topo- 
logical group T" whose elements label the distinct states. In 
our models, the definition of this loop operator is trivial and 
transparent: it is just the generalized "height difference". 



B. Outline of the paper 

In this paper, I first (Sec. [TTJ) generalize the height-model 
idea to the case where the height variable belongs to a discrete 
(finite) abelian or non-abelian group, thus defining a family 
of classical models which (in many cases) has a topological 
order. The models are defined by a lattice, a group, and the 
selected subset of group elements which are permitted values 
for the 'spins" of the model; the spins sit on the bonds. Non- 
Abelianness of the group has interesting consequences: for 
example, a collection of defects no longer has a unique net 
charge (Sec.lITCl) 

In Sec. [HI] and Sec. IIVI I survey the various combinations 
for the smallest non- Abelian groups, using the crude Pauling 
approximation as a figure of merit to identify the most attrac- 
tive models for Monte Carlo simulation, using single-site up- 
dates. Furthermore, Sec.[V]suggests what quantities are inter- 
esting to measure in such a simulation; however, no simula- 
tion results are reported in this paper. But some first analytic 
results are included in Sec. lVIl based on transfer matrices and 
hence implicitly one dimensional: the main point is to shed 
light on how the size dependence or defect pair correlation 
depends on the group elements labeling the topological sector 
or the defects. 

Finally, the conclusion (Sec lVHt reflects on which topolog- 
ical behaviors are inherently quantum mechanical, and which 
are not (in that the same behavior can be found in classical 
models). Furthermore, applications are suggested, either to 
simulating systems with vacancy disorder, or to constructing 
quantum versions of the models in this family. 

n. DEFINITIONS AND TOPOLOGICAL BEHAVIORS 

This paper is meant to introduce (and compare) a whole 
family of models. In this section, I define the general rules for 
this family, and then describe the most promising examples. 
(In the next section, I shall exhibit consequences for Monte 
Carlo updating of such models.) 

A. Model definition 

Let us take a "lattice" (not necessarily Bravais, e.g. honey- 
comb) of sites r. The spins sit on the bonds of this lattice, and 
take values in the discrete group Q ; they are 

<x(r,r')e£ (2.1) 

where (r, r') labels a bond of nearest neighbors. The bonds 
are directed; if we reverse the direction we view a bond, the 
spin on it turns into its inverse: 

cr(r',r)=a(r,r')-\ (2.2) 

Then each configuration of the spins induces a configuration 
of "heights" h(r) E Q, defined by 

h(r) =a(r,r')*h(r'), (2.3) 
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where "*" represents the group multiplication. Of course, 
h(r) is only defined modulo a global multiplication by some 
element r, h'(r) — h(r)r; to make it well-defined, we could 
arbitarily require (say) h(0) = e, where e is the group iden- 
tity. One could then explicitly construct h(r) at the neighbors 
of site 0, and iteratively at their neighbors, etc.; the result is 
independent of which bonds (r,r') are used for this, if and 
only if a plaquette constraint is satisfied [Eq. i2.5i . below]. 

It will be useful throughout to define the line or loop prod- 
uct of p spins: 

7(f) ee cr(ri,r p ) * a(r v ,T p -i) * ... * cr(r 2 ,ri) (2.4) 

Here the loop £ is a string of bonds (rfc, r^+i) connecting end- 
to-end, for k = 0, ...p; it is a loop when r p = tq. 

1. Plaquette constraint 

Two constraints are imposed on the spin configurations. 
The first is the plaquette constraint: we require the loop prod- 
uct around any elementary plaquette, to be the identity, 

7(Vq) = 6 (2-5) 

This is necessary (and sufficient) for h(r) to be well-defined. 

One can define variants of any model by relaxing the pla- 
quette constraint to allow a small number of defect plaquettes, 
around which the loop product is not the identity. Suppose 
that defects cost an energy A (possibly depending on the kind 
of defect): then the basic, defect-free version of the model can 
be viewed as a Boltzmann ensemble in the limit T/A — » 0. 
On the other hand, if we imagine there were a spin-spin inter- 
action that breaks the degeneracy of the states satisfying the 
plaquette constraint, then the basic version of the model is the 
T/J oo limit. 

Using condition (12. 5t and induction (adding one plaquette 
at a time to the loop), the loop product must be j(£) = e for 
any finite loop £ that is contractible to a point in small steps. 

2. Spin constraint 

The second constraint is the spin constraint: choose a "spin 
subset" S C Q such that 

a(r,r')eS (2.6) 

everywhere. Hence, the choice of S is a major part of a 
model's definition; in Sec. IIIICI below, I will discuss other 
desirable features of S. The spin constraint is retained even in 
versions of the model with defect plaquettes. 

The constraints should respect the group and lattice sym- 
metries . Implementing the group symmetry means requiring 

(re5=^u*cr*ii _1 e S (2.7) 

for any conjugating element u. Thus, S must be one of the 
group's conjugacy classes - the simplest case - or a union 



of such classes. (Some non-abelian groups, and all abelian 
ones, have "outer" automorphisms, symmetries which cannot 
implemented by conjugacy within Q: we may also wish to 
implement those symmetries, too). 

To implement lattice symmetry, one asks 

a e S =► a' 1 e S (2.8) 

so the model respects inversion (around the bond's midpoint). 
H 

Let us further require 

e i S (2.9) 

thus a uniform height configuration is disallowed. Finally, and 
trivially, 

S generates the full group Q. (2.10) 

(If not, I could have redefined Q as the subgroup generated by 
S.) 

Without the spin constraint, the models would be identical 
to the lattice gauge models of Doucot and Ioffe (23J. In such a 
model, only gauge-invariant (i.e. loop) quantities can have 
nonzero expectations; other correlation functions are zero, 
even at the nearest distance. In contrast, these group-height 
models (like the original height models) have nontrivial finite- 
size effects and local correlations: in particular, there are me- 
diated interactions between topological defects. 

Furthermore, the spin constraint allows the possibility of a 
long-range ordered phase, particularly if we assign different 
Boltzmann weights to different configurations, and we may 
find phase transitions as those parameters aare varied. Partial 
orderings are also possible, and transitions might occur be- 
tween different topological orders It will be easier to explore 
the phenomenology of such transitions in the classical realm. 

B. Topological sectors and topological order 

For these models, a "sector" means simply the configura- 
tions that can be accessed by a succession of local updates. 
(Here "local update" means an operation that turns one valid 
configuration to another by changing spins in a small neigh- 
borhood of some site, as might be deployed for Monte Carlo 
simulation. A "nonlocal" rearrangement, as developed in 
Sec. [IV] means the cluster of updated sites can be arbitrarily 
large, and in particular could include a topologically nontriv- 
ial chain of sites that spans the periodic boundary conditions.) 
By this definition, "sectors" are well defined in any finite sys- 
tem larger than the maximum update cluster. In other models 
and with other definitions, passing between sectors is be abso- 
lutely forbidden, so that sectors (more exactly "components", 
like the up and down ordered phases in ferromagnet) are emer- 
gent in the thermodynamic limit. 

In these models, sectors can be labeled by loop products 
7(£). As noted after J2.51 l, products around topologically triv- 
ial loops must give the identity, but others - e.g. through the 
periodic boundary conditions of a torus system - in general do 
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not. Such loop products are not changed by local updates, and 
therefore must take the same value for all states in a sector. 
So we call sectors "topological" when they are distinguished 
(and necessarily disconnected) by having different values of 
the loop product(s). The definition of topological order is that 
- in the thermodynamic limit - different topological sectors 
all become equivalent, in that they cannot be distinguished by 
any local expectations. Furthermore, just as the ground state 
energies in different sectors should become equal in the case 
of quantum topological order, the free energies should become 
equal in our models. 

Let £1,(2, ■■■(■g be the basic independent loops, where g is 
the genus; then the loop products {7^ = 7^)} [all taken 
from the same origin] label the possible topological sectors 
of configuration space. An interesting question is how many 
distinct sectors there are, [24] given the system's genus g. 

1. Invariance 

Before counting sectors, we need to explore invariance 
properties of the sector labels, in case some labels are equiv- 
alent to others. First, there is a sort of gauge freedom: if we 
had evaluated these loops starting from r instead of from 0, 
then 

It -> it = 70r * le * lor 1 , (2-11) 

where 70 is the line product along any path from the origin 
to r; notice that this same element conjugates all the distinct 
loops. However, just because our labeling fails to distinguish 
two sectors does not conclusively show they are the same. 

A better criterion for counting sectors as equivalent is 
that one can be turned into another by local updates. Keep 
the same origin, but perform a single-site update [see 
Eq. eqreq:sigma-new-inner, below] hitting on the origin ver- 
tex, one gets another conjugacy 

It -> it = t*7£*t _1 , (2.12) 

where t is now the updating multiplier. 

When the group Q is abelian the sector labels are invariant 
with respect to how we take the loop and unchanged by local 
updates. We can have an independent and invariant loop prod- 
uct 7i for every topologically independent loop, so the number 
of sectors is ng 2g where g is the system's genus (2g = 2 for 
torus), and ng is the number of elements in the group. 

2. Sector counting in non-abelian case 

On the other hand, in the non-abelian case, a loop product 
is invariant only up to a conjugacy, so we have fewer sectors. 
Furthermore, the allowed values of distinct loop products are 
not independent. Consider a square lattice model in a rect- 
angular system cell L x x L y , with periodic boundary condi- 
tions; let (j x , j y ) be the loop products along straight lines of 
bonds running from the origin site (0, 0), in the x or y direc- 
tions respectively. The loop running from (0, 0) to (L x , 0) to 
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FIG. 1 : The composite defect charge of a pair may be changed by 
sliding a third defect between the two. (a). Two defects (stars) have 
charges 71 and 72, thus a loop enclosing both contains charge 7271. 
A third defect with charge 73 is being moved from the right. The loop 
product around defects 3 and 2 is 7273 These charges are defined 
using the reference point P, and multiplications are from the right, 
(b). After defect 3 is moved between defects 1 and 2, the product 
around defects 2 and 3 is still 7273. (c). After defect 3 has passed 
to the other side, the product around defect 2 has the new value 72; 
the product around 2 and 3 is 7372, which is unchanged from (b) and 
therefore equal to 7273. Hence 72 = 73 72 7^ > an d the combined 
charge of defects 1 and 2 is now 7271 = 73727 3 _1 7i, which can be 
in a different class than 7271. 

(L x , L y ) to (0, L y ) and back to (0, 0) contains no defect, so 
by inductive use of the plaquette contraint its loop product is e. 
Yet the four segments of this loop are just j x and j y , forwards 
or backwards, so the loop product is 

ly 1 *"fx 1 *Jy *1x = e; (2.13) 

in other words, j x and "f y must commute. 

So, in effect, we must define an equivalence relation 
(lx,1y) ~ (ixiiy) whenever the pair satisfies ( |2.13t , and 
each topological sector is one equivalence class. In the case 
of a larger genus, we extend in the obvious way to longer lists; 

Some obvious kinds of equivalence classes are: 

(i) (e,e) 

(ii) (u,e) or (e,w) 

(iii) (oj, oj k ) or (cu k , uS) for k = 1, m — 1, where ui is an 
element (not the identity) of order m. 

(iv) If the group has a nontrivial center Qz, consisting of 
elements that commute with all the other elements, then if 
{lx! ly) is a sector then (z x -f x , z y -f y ) is another sector, where 
Zx,z y S Qz- 

TableQ]shows the number of classes /Ui and the sector count 
fi2 for some groups of interest. 

C. Defects and non-abelian effects 

We can allow the possibility (dilutely) of a plaquette that 
violates the plaquette constraint. The loop product around it 
will be called f3. It is analogous to the Burgers vector of a 
dislocation, or of the topological defects in the usual height 
models based on Z). 

In the case of a height model (in a rough phase), topolog- 
ical defects are like vortices in a two-dimensional Coulomb 
^as lT2lll22Tl . They behave like U(l) electric charges in a two- 
dimensional universe. Opposite charges feel an attractive log- 
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arithmic potential. In contrast, in the case of topological order, 
the attractive potential decays exponentially and the defects 
are deconfined B27I1 . 

If the topological order is non-abelian, the non-commuting 
property of defect charges has some interesting consequences. 
They are not unique to topological order; this also is a long 
known property of defects of traditional ordered states asso- 
ciated with a non-abelian homotopy group 12811 . One conse- 
quence is that a loop product may be changed when a defect 
of charge j3 is passed across it, the action being a conjugation: 

7(£)^/?7W/r 1 . (2.14) 

Thus the topological sector might be changed when a defect 
wanders around the periodic boundary conditions. Also, the 
net charge of a defect pair can be changed by passing another 
defect between the pair. (See Fig. [TJ 

Another consequence of non-abelianness is that given two 
given defects of specified charges, there is more than one 
possible value for their combined charge. In the quantum- 
mechanical approaches to defects in topologically ordered 
systems, this same property is also the hallmark of non- 
abelianness. In that context, the list of allowed combinations 
is known as the "fusion rules", and there are matrices (gener- 
alizations of Clebsch-Gordan coefficients) which tell how to 
form the appropriate linear combinations. 

A third consequence is pertinent to simulations and the def- 
inition of topological sectors in the presence of defects. If 
one creates a defect pair and moves one defect around the 
boundary conditions, it may recombine with the original de- 
fect into a single defect, rather than annihilate. The fact that 
two defects may not be able to re-annihilate is very similar to 
the "blocking" idea of Ref. ll24ll (for quasiparticles in a non- 
Abelian quantum Hall state). 

The single defect state satisfies the generalization of ( 12.131 ), 
namely 

ly 1 * Ix 1 * ly * Ix = P- (2.15) 

This commutation is a "group commutator". Such a single- 
defect state may conveniently allow numerical measurements 
of the creation free energy of a single defect. Of course, such 
a state is never possible for abelian defects; in that case, a 
system with periodic boundary conditions must have either 
no defects, or at least two of them. 

III. POSSIBLE MODELS 

In this section, I survey specific models, emphasizing the 
criteria which would make some of them particularly attrac- 
tive for future investigations. To summarize Sec. Ill Al models 
in this paper are specified by (i) the group Q (values of height) 
(ii) the spin subset S (values of spins) (iii) the lattice whose 
bonds the spins sit on. 

Therefore, the models will be named in the form 
"<?(m)latt". Here "G" is the groups name, (m) is the or- 
der of the elements in the selected conjugacy class (usually 
that is unambiguous), and "latt" abbreviates the lattice ("tri", 



product y x = b 




(a) (b) 

FIG. 2: [Color online.] (a) Example configuration of the (abelian 
group) model Zi x Z2(2)sq. Each square lattice edge is occupied 
by a group element a, b, or c; directions are unneeded since each of 
these is its own inverse. The loop products and 7 y around the pe- 
riodic boundary conditions are also shown, which define the topolog- 
ical sector. For an abelian group, their values are independent of the 
starting points, (b) Example configuration of the (non-abelian group) 
model 53(2, 3)tri. Labels a, b, c denote the elements (23), (13), (12), 
while the arrow denotes a cyclic exchange (132). 

"sq", or "he" for triangular, square, and honeycomb). Thus 
"5*3(2, 3)tri" means that Q is the permutations of three ob- 
jects, S contains all three pair exchanges, as well as the two 
cyclic permutations (i.e. every group element except for e), 
and "tri" means the spins sit on the edges of the triangular 
lattice. An variant nomenclature is sometimes convenient, in 
which the "(m)" in the label gets replaced by "{<7i, a%, ...}": 
the set {<7i,<72, ■■■} is simply the listing of the selected ele- 
ments. 



A. Groups 

Table|I]lists the groups and spin subsets I shall be interested 

in. 

For future reference, I mention the automorphism group 
Aq of a group Q, which is simply its symmetry group. Each 
a e Aq is a permutation of the group elements preserving its 
structure, a(gg') = a(g)a(g'). When Q is non-abelian, there 
is a subset of the automorphism group called the inner auto- 
morphisms, defined as the conjugations, a T (g) = rgr^ 1 . Ob- 
viously the inner automorphism subgroup 
is isomorphic to G/Gz, where Gz (the center subgroup) con- 
sists of the elements that commute with everything. But many 
groups have additional outer automorphisms that are not con- 
jugations; in particular, all automorphisms of an abelian group 
are outer. 



1. Abelian groups 

We start by considering discrete abelian groups in this fam- 
ily of models. The smallest of them, Z2, does not work since 
S can have only one element. The next simplest cases are 
cyclic groups Z q , i.e. the integers modulo q under addi- 
tion, although these often turn out to be height models (see 
Sec.|IIIB2Jbelow). 

Beyond that we go to direct products of cyclic groups, in- 
deed any abelian group can be represented thus. If S was also 
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TABLE I: Groups and spin subsets, /i is the number of conjugacy 
classes, and p,2 the number of topological sectors on a torus; ng, ns, 
and ne respectively are the number of elements in the group Q, the 
number in the selected subset S, and the number of even elements. 
The effective bond probability pb is given by formula J4.4t . 



group+tag 


ng 


/' 


M2 


n s 


ns 


Pb 


Z 2 x Z 2 (2) 


4 


2 


16 


3 




1/3 


S3 (2) 


6 


3 


8 


3 


3 





S3 (2,3) 








5 




1/5 


Q(2) 


8 


3 


10 


6 




2/7 


D4{m, m'} 


8 


5 


20 


4 


4 


4/9? 


AtO) 


12 


3 


8 


8 




4/11 


A s (2) 


60 


4 


20 


15 




45/59 


As(3) 








20 




40/59 


A 5 (5) 








12 




48/59 



taken to be a direct product, of course the model would reduce 
to a superposition of non-interacting models, one for each fac- 
tor. However, there are many attractive examples in which S 
is not a direct product, in particular Z 2 x Z 2 (see Sec. lIIIBTt . 



2. Non-abelian groups 

The smallest non-abelian group is S3, the permutations 
on three objects (also isomorphic to the dihedral group D3). 
Here, S may be taken as the class of all pair permutations, the 
model 6*3(2), or as all permutations except the identity, that is 
63(2, 3) in our notation. 

Each of the next two smallest non-abelian groups has eight 
elements. One of these is the 8-element quaternion group Q, 
i.e. the unit elements {±1, ±i, ±j, ±k} from the quaternion 
ring. Here S must be the class of the six elements not equal to 
±1. 

The other eight-element non-abelian group is D 4 , the sym- 
metry group of the square lattice. 

The "alternating groups" A4 and A§ are especially attrac- 
tive for our purposes due to their high symmetry (so we can 
choose S to be a single class containing a sizeable fraction 
of all the group elements). They consist of the even permuta- 
tions of four and five elements. Note that A4 and A$ are also 
the point groups of the (proper) rotations of a regular tetra- 
hedron and a regular icosahedron, respectively. Being sub- 
groups of SO(3), these groups might in some sense serve as a 
discretization of it 12911 . just as clock models are a discretiza- 
tion of the XY model. That would be interesting as a way to 
make a connection to topological models (or gauge theories) 
defined in terms of Lie (i.e. continuous) groups. 

Finally, A$ is the smallest non-abelian simple group, mean- 
ing it has no normal subgroups; as we shall see in a moment 
(Sec. IIH Bl i. normal subgroups are an annoyance since they 
tend to make the behavior more trivial than would be expected 
for the group Q. 



B. Example models 

Next I shall survey the simplest examples. Most of them 
reduce, in some fashion, to previously known models; that is 
an advantage for computational studies, since old results can 
be used as checks. In several cases, the models in our family 
have "accidental" topological order, i.e. beyond the group Q; 
In particular, some of them have height representations. 

The group and subgroup involved in our spin constraint are 
finite, and so is each plaquette; thus it can happen that the 
allowed configurations satisfy stronger constraints than those 
they were designed to fulfill. The first five subheadings below 
all, in one sense or another, reduce to known models. 



1. Z2 X Z2 and the 3-coloring model 

For a first example, let Q = Z2 x Z2, an abelian group. 
Besides the identity, this group has three equivalent elements 
a, b, c; each has order two, and the product of any two gives 
the third. If we treat these as a class (although they are not 
conjugate, since the group is abelian), then we must choose 
that class to be the spins, Z 2 x Z 2 (2). Since a = a -1 , etc., we 
can depict the spins using three (undirected) "colors" of the 
edges. On the square lattice this gives perhaps our simplest 
example (Figure |2). 

What about the triangular lattice case [model Z 2 x 
Z 2 (2)tri\l The plaquette constraint is simply that each trian- 
gle has one edge of each color: the "three-coloring model". (It 
is usually represented on the edges of the dual [honeycomb], 
where the constraint says each vertex has three colors; in ei- 
ther case, the spins live on kagome lattice vertices, and the 
configurations are also the ground states of the 3-state Potts 
antiferromagnet on that lattice.) This model is known to have 
aZxZ height representation IU7ll26ll . in addition to the finite- 
group height field h(r) defined by ( 12.31 ). 



2. 6-vertex model 

For another example, take Q to be Z q , with q > 4, and let 
the lattice {r} be the square lattice. Choose S = { + 1,-1}. 
(The two elements are not the same class; they are related only 
by an outer automorphism.) Then the sum of spins around a 
plaquette can be zero (mod q) only if it is just zero, i.e. there 
are exactly two +1 and two —1 in the loop. If we express these 
spins on the dual (also square) lattice, as an arrow pointing 
outwards (resp. inwards) wherever u = +1 (resp. —1) as 
the loop is traversed counterclockwise, we see these just are 
the configurations of the six-vertex model - which also has a 
integer-valued height field. Ir42ll 

Since Z or Z m are abelian groups, {+1, —1} is merely an 
outer class. 
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3. Groups with an even subgroup 

An even subgroup £ (with ri£ = rig/ 2) has Q/£ = Z2. 
That is, any product of an even number of elements lies in £ . 
Say that the spin subset S consists of odd elements. (If it con- 
sisted of even elements, we would generate at most £.) Notice 
that such a model cannot use the triangular lattice, since the 
plaquette rule cannot be satisfied (the plaquette product must 
be odd, while e is an even element). 

Now if the simulation cell has even dimensions, the possi- 
ble topological products 7(4) must lie in £. (Even if the cell 
has an odd dimension, the possible values of j(£i) still corre- 
spond 1-to-l with elements of £ .) Thus, the topological sector 
labels can only belong to the subgroup £. 

For example, in permutation groups the even subgroup con- 
sists of even permutations. In the case of S3, the even permu- 
tations are just (123) and its powers, so £ = Z3. Conse- 
quently the model 5s(2)sq can have only abelian topological 
behavior. In our list, D4 is another group that contains an even 
subgroup (the proper rotations). 

4. Groups with a center 

What if Q has a non-trivial center Qz 1 - (The center is sub- 
group of elements that commute with every other element). 
For example, if we adopt the group Q of unit axis quater- 
nions, (which order 8) then Qz = { + 1,-1} = Z2 and 
Q/Qz — Z2 x Z2- Thus, the model Q(2)tri projects onto 
configurations of the the three-coloring model. (Just map 
±i — > a, ±j — > b, ±k — > c.) The group D± also has a cen- 
ter (two-fold rotations, i.e. inversions, commute with every- 
thing.) 

C. Criteria for models: estimates of entropy 

To estimate at once the viability of many different models, I 
shall use very crude estimates of the entropy and (in Sec. lIV B"b 
updatability. Say the lattice has coordination z and the dual 
lattice has coordination Zd, i-e. the number of sides of each 
plaquette. (These numbers are related by 1/z + l/zd = 1/2.) 
Also, say the group has ng elements, of which ng are in the 
selected subset S. These three parameters — ng, ris, and z 
(or Zd) - contain much of what we need to characterize the 
possible models. See Table iHl for the parameters related to 
lattice geometry, and Table Q] for those related to the groups 
and the spin subsets. 

I will use a Pauling estimate for the entropy. There are 
Nz/2 edges and hence ng Nz ^ 2 ways of placing spins inde- 
pendently chosen from S, in a hypothetical ensemble that does 
not (yet) enforce the plaquette constraint. If we knew the frac- 
tion of all these states that do obey the plaquette constraint, we 
would have the total count of allowed states and thus the de- 
sired entropy. 

Pauling's approximation is to pretend the event of satisfy- 
ing the plaquette constraint is uncorrected between plaque- 
ttes. So let / e be the chance that a given plaquette has plaque- 



TABLE II: Lattices (asterisk denotes an average over two kinds of 
plaquettes). Here "rr-phase lattice" denotes the lattice (3 2 ,4, 3,4) 
and "square-octagon" lattice denotes (4, 8 ). The columns give the 
coordination numbers z of the lattice and Zd of the dual lattice, fol- 
lowed by the lattice's bond and site percolation thresholds p c b and 
p cs . (Many more significant digits are known 13411 .) 



lattice tag 


Z Z d Pcb Pes 


triangular tri 
cr-phase 

square sq 
kagome 

honeycomb he 
square-octagon - 


6 3 0.347 0.500 
5 3.333* 0.414 0.551 
4 4 0.500 0.593 
4 4* 0.524 0.653 
3 6 0.653 0.697 
3 6* 0.677 0.730 



tte product equal to e. Then in this approximation, the prob- 
ability is j^ z l Zd that the plaquette constraint is satisfied on 
all Nz/zd plaquettes of the whole system. Thus the Pauling 
estimate of the ensemble entropy is 

e iVS PauIing _ ( s z/2jz/z d \ (3 1} 

The condition we must satisfy, in order to have an ensemble 
at all, is Spauiing > 0, i.e. 

n s > l// e 2/2d ., (3.2) 

If Zd is not too small, we may estimate that the plaquette 
product is equally likely to be any group element, hence very 
crudely 

fe « l/ng. (3.3) 

Less crudely, one can work out the the actual probabilities that 
the product of Zd random group elements from the allowed set 
S will give the identity, and these are the f e values in TablelHll 
Comparison of those f e values with ng in Table Q] shows that 
usually, d3.3t is not bad. When the product of two elements 
of S is particularly likely to fall into one class, the true f e de- 
viates more from d3.3t , either on the low side (e.g. the model 
A5(2)tri) or the high side (e.g. As(5)sq). The extreme case is 
if the group contains even and odd elements, and S consists of 
odd elements (the 6*3(2) or D±(m, to') examples in TablelTTl]). 
In that case we should replace ( 13.3b by f e rj 1 /ng = 2 /nG if 
Zd is even, but f e = if Zd is odd. 

There is just one entry in Table U showing a negative Paul- 
ing entropy Spauiing < 0, namely A5(2)tri. It is convenient to 
explain this case in the language of (proper) rotation group of 
an icosahedron, which is isomorphic to A5. The only way that 
three twofold elements can multiply to give the identity is mu- 
tually when the two fold axes are mutually orthogonal. Since 
the triangles share edges, any valid global configuration must 
use that same triad in every triangle; this entails a fivefold 
(Sg) global symmetry breaking, since the fifteen twofold axes 
of icosahedral symmetry break up into five disjoint orthogo- 
nal triads. Indeed the three used elements form a subgroup 
isomorphic to Z2 x Z2 so we are back at the three-coloring 
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model for which <Sp a uiing > 0. The Pauling approximation 
gave zero entropy only because it did not take account of the 
symmetry-breaking and attempted to mix domains with in- 
compatible symmetry breakings. 

My purpose here is not to obtain quantitative estimates of 
the model's entropy, although the Pauling estimate is some- 
times surprisingly accurate. Rather, I want to compare these 
values between different models as a figure of merit, to aid 
us in guessing which models are the most interesting or the 
most tractable. To this end, the figures of merit are shown in 
Table Unl 

To satisfy Eq. 13.21 the three parameters get pushed in the 
following directions, but there are considerations limiting 
each of the three. 

(1) We want ng as small as possible; however, there are 
not so many small, discrete, non-abelian groups: only 
three have ng < 8, namely S3 (permutations of three 
objects), Q (quaternion group), or D4 (point group of a 
square). 

(2) We want larger ns, meaning the model is less 
constrained (and more tractable). In the limiting case 
ns = ng, the model is just a pure gauge theory, which 
is trivial apart from its global topological properties. 
On the other hand, a sufficiently large ns requires in- 
cluding more than one conjugacy class in S, so that the 
spins can have inequivalent "flavors". That is estheti- 
cally undesirable: a generic model (with unequal statis- 
tical weights) needs more parameters, and it is harder to 
imagine how such a model could be realized physically. 

(3) We want large Zd, as in the honeycomb lattice. 
However, it is esthetically harder to implement a prod- 
uct constraint in a physical model. (When the prod- 
uct string is short, there are only a few symmetry- 
inequivalent cases for it, and it is easier to concoct a 
Hamiltonian term which does not reference the group 
multiplication, but which has those cases as its energy 
minimum.) 



To satisfy (13.21 with a large group but S consisting of just 
one conjugacy class, the group must have high symmetry. 
E.g., the alternating group A$ (the proper icosahedral rota- 
tions) has ng = 60 and contains conjugacy classes with 12, 
15, or 20 elements, which using (2) would need Zd > 3.30, 
3.02, or 2.41 respectively. 



IV. MONTE CARLO UPDATING 

For us, one essential criterion of a model is the possibil- 
ity of Monte Carlo simulation. I limit consideration to the 
equal-weighted ensemble, in which every allowed configura- 
tion has the same weight. Then detailed balance is satisfied 
if the forwards and backwards rate constants are the same for 
any update move. But what is the minimum sufficient update 
move? For the six-vertex model it sufficed to reverse the ar- 
rows on the four edges of one plaquette, which changes the 



height field on one site, a purely local update. For the three- 
coloring model, the minimal update involves switching two 
"colors" (e.g. a <H> b) along a loop, a nonlocal update move. 
What happens generic ally for our family of models? 



A. Cluster update move 

The update move is simplest described in terms of the 
height function (defined in ( 12.3b ). First pick at random a group 
element T^e and a starting site ro. Say T> is the domain be- 
ing touched by the update. (It will be explained in a moment 
what determines T>). Then I prescribe that the update premul- 
tiplies the heights in this domain by r, so as to "shift" them: 



h'(r) = 



t * /i(r) for reP; 
h(r) for r (£ V. 



(4.1) 



This induces the following update of the spin configura- 
tion: H 



a'(r',r) = 



T*a(r',r))*r 1 forr,r'eP; 
r*cr(r',r)) forr' eX>,r' £V; 

o-(r',r)) * r" 1 forr' (£"D,r' eV; 
cr(r', r) for r, r' V. 

(4.2) 

I call this a "gaugelike" transformation 13311 : it has the same 
form as a gauge transformation would, but it is valid only 
when an additional spin constraint is satisfied too. 

If both endpoints of the bond are in T>, then a' is conjugate 
to a and must be legal (since we include whole conjugacy 
classes in S). 

On the other hand, where the (r,r') bond crosses the do- 
main boundary &D, the spin constraint is nontrivial to satisfy. 
Let's place an arrow along the edge from r to r' if and only if 



(4.3) 



In other words, there is an arrow from r to r' whenever in- 
cluding r in T> forces us to include r' as well. This arrow is 
not bidirectional, (it is in the case t 2 = e). Thus, we might 
have any of four possibilities (no arrows, arrows both way, or 
arrows one way) along each bond. 

Then the update rule is to construct the arrowed-percolation 
cluster consisting of site ro, with the rule that site r' is in- 
cluded if site r is included and there is an arrow r to r'. This 
is the smallest possible updated domain containing ro. Of 
course, we do not actually need to construct all the arrows; in- 
stead, we grow the cluster from the initial site, and construct 
arrows only from sites already in the cluster. 

Notice that (only) in the case Q is abelian, ( 14.21 ) reduces 
to cr' = a throughout the interior of T>. In other words, the 
update only changes spins along the boundary dT> and thus is 
a loop update. In a non-abelian model, however, the update is 
generally a cluster update. 

In some models (see next subsection) there is a strong 
chance to hit a system-spanning cluster, including most of the 
sites, which tends to be inefficient. (Updating all the sites is 
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equivalent to no update). To avoid this, a limiting size s max 
for the update cluster T> should be set; if this limit is reached, 
we cancel the tentative move and start over, choosing a new 
random ro and r. 

B. Numerical criteria for cluster updates 

Notice that in growing a cluster from r , we never cared 
about the reverse arrows. Therefore, we obtain the same clus- 
ters as we would in an ordinary (not arrowed) percolation 
problem, if the occupied bond probability pb is identified with 
the probability of an arrow in a pre-selected direction; that 
probability is simply 



if we chose the candidate updating factor t at random. These 
probabilities are shown in Table Hill In a case that any group 
element r works as a multiplier on any bond, I would write 
Pb = in Table HIH rather than use (14.41 ). In such a case, our 
model is locally trivial: exactly ng N configurations may be 
accessed, simply by applying one arbitrary group element at 
every site. In other words, there is a (locally) 1-to-l mapping 
to the trivial model in which every site has an independent 
degree of freedom. That model is just the gauge model, which 
was studied previously in Ref. [23]. 

In the spirit of the Pauling approximation, let us now pre- 
tend that arrows on different bonds are uncorrelated: within 
that assumption, we must obtain the same cluster distribution 
as in the (thoroughly studied) problem of uncorrelated perco- 
lation on these lattices. It follows that the updating behavior 
tends to depend on the relation of pb to the critical percolation 
fraction p c b- On the one hand, if pb > p c b, then the cluster 
grows without limit, including a nonzero fraction of the whole 
system; in that case, the update move certainly is not efficient. 
On the other hand, if pb /p c b is too small, we never get a clus- 
ter at all, or else a single-site update (next subsection) would 
suffice. The "interesting" case when a cluster update is neces- 
sary and helpful, would be for Pb/p c b close to or slightly less 
than unity. 

C. Single-site updates 

A single-site update is the case that the updated cluster T> 
is just one site, thus only the z spins around it are updated. 
When pb is much less than p c b, most clusters are small, and 
the probability Pi of a single-site update is appreciable. If 
Pi is large enough, it might be ergodic to use only single-site 
updates (i.e. to pick s max = 1), in which case we can omit 
the cluster-growing algorithm. I shall concentrate on these 
cases, which are the easiest to simulate (and also the likeliest 
to extend to quantum models). 

To estimate Pi, I pretend the z bonds around a site are inde- 
pendently occupied by randomly chosen elements of S. Then 

Pi est = 1- JJ(1- (4.5) 

a 



TABLE III: Entropy and updatability parameter estimates for se- 
lected models. Formulas from eqs. d3.ll >. J4.4I ). and ( 14.51 ). Note a: in 
these cases, any even element r can always update, but no odd r can 
ever update. 



Model name 


fe exp(S 


Pauling) Pb/Pcb 


Pi est 


Z 2 x Z 2 (2)tri 


2/9 


4/3 ... 


0.241 


Z 2 x Z 2 (2)sq 


1/3 


7/3 0.667 


0.681 


S 3 (2)sq 


1/3 


3 0.0 


0.5 a 


S 3 (2)hc 


1/3 


3 0.0 


0.5 a 


S 3 (2,3)tri 


4/25 


16/5 0.576 


0.870 


S 3 (2,3)sq 


21/125 


21/5 0.400 


0.963 


0(2)sq 


7/54 


14/3 0.571 


0.930 


D4{m, m'}sq 


1/4 


4 0.0 


0.5 a 


D4{m, m'}hc 


1/2 


4\/2 0.0 


0.5 a 


A 4 (3)tri 


1/16 


2 0.727 


1.000 


A 4 (3)sq 


3/32 


6 0.613 


1.000 


A(2)tri 


2/225 


4/15 2.20 


0.668 


4 5 (3)tri 


7/400 


49/20 1.95 


0.894 


A 5 (5)tri 


5/144 


25/12 2.34 


0.706 


A 5 (2)sq 


71/3375 


19/5 1.53 


0.285 


A 5 (3)sq 


147/8000 


147/20 1.36 


0.544 


yl 5 (5)sq 


53/1728 


53/12 1.63 


0.360 



where a indexes the group class (with n a elements) that r 
might be in (excluding the identity), and I defined q a to be the 
fraction of times t * a E S given that r falls in class a. 

To digest the implications of ( 14.51 ). let's make an even 
cruder version of the estimate, replacing q a by it's average 
over all r's, namely q a —> 1 — pb'. I get 1 — [1 — (1 — Pb) 2 ]™ 5-1 
which is a lower bound on Pj 08 ' as given by ( 14.51 ). Evidently, 
to have a high single-site success rate, we want (i) ng as large 
as possible, (ii) z as small as possible, and (iii) pb as small 
as possible; in light of ( 14.41 ). the third criterion amounts to 
wanting n$ /ng as large as possible. Those are the same three 
considerations given in Sec. Hild as favoring a large entropy. 

I include these estimates in Table [Till particularly focus- 
ing on the models using group A 5 . We see from Table Hill 
that p^ 8 * is large enough in many cases that we can rely on 
single-site updates. However, whenever p^ gets close to 
1, our model is "too easy" in some sense - it is practically a 
gauge model, with only mild constraints eliminating some of 
the configurations. 

The entry Pi° st = 1 for At (3) is delusory. This comes be- 
cause q a = 1 for a certain class of update multipliers, namely 
the order-2 class (double pairwise exchanges). If we limited 
ourselves to this class, indeed every update would be success- 
ful, but (it can be checked) the move would not be ergodic 
(does not access the whole ensemble). A similar situation ap- 
plies in the cases of £3(2) or D 4 {m, m'}: any r from the 
even subgroup is always accepted, while an odd r is never ac- 
cepted; but single-site updates based on the even subgroup do 
not access the whole ensemble. 

To implement an actual simulation, one would not want 
to choose r at random, but biased towards the group classes 
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with a larger q a (the success fraction looking at just an iso- 
lated bond). In particular, one would omit group classes with 
q a = 0; if the group contains even/odd elements and S in- 
cludes only one parity of element, then q a = for every class 
of odd elements. The values of pt, and p^ 81 in Table UTTIfor 
53(2) and £>4(m, ml) were computed assuming r € £ . 

Another way to implement a single-site update is, after 
choosing a random vertex r, to examine the local environment 
of its z bonds, find the entire list of r's which can update it, 
and choose randomly from this list. Typically, configuration 
dependent choices like this are avoided in Monte Carlo algo- 
rithms because they tend to violate detailed balance. In the 
present case, however, it can be checked that the number of 
possible t's is always the same in the old and new configu- 
ration, i.e. the rate is the same for the forward and backward 
step, which is sufficient to ensure detailed balance (and an 
equal ensemble weight for every configuration). 



D. Criteria for initial conditions 

In height models, certain special states (e.g. the "columnar" 
arrangement of dimers on the square lattice) were "ideal" in 
having the maximum number of possible update moves. (For 
a model requiring loop updates, we might replace that cri- 
terion by "having the shortest typical loops.") Certain other 
states (e.g. the "herringbone" packing of dimers) were inert, 
in that no finite updates are possible (in the thermodynamic 
limit). These states, in a height model, correspond respec- 
tively to a zero coarse-grained gradient of the height variable, 
or the maximum gradient. 

In a non-abelian height model, the coarse-grained height 
gradient is undefined, but one can still construct "ideal" and 
"anti ideal" states. It is recommended that simulation runs be 
started in both kinds of state, being in some sense opposite 
extremes of the configuration space. A diagnostic for equi- 
libration is then whether the expectations from the two starts 
converge to the same values. 

More exactly, rather than a single domain of anti-ideal state, 
one should divide the system into two domains. Then, updates 
are initially possible along the domains' border, using loops 
which extend across the system. Gradually, a larger fraction of 
the system's area become updatable, and the loops get smaller. 
On the other hand, starting from an ideal state, the loops are 
initially small and get larger. Thus, tracking the loop distri- 
bution is an obvious diagnostic to test for convergence to the 
same equilibrium state. 



V. POSSIBLE MEASUREMENTS IN SIMULATIONS 

In this section, I sketch how one might confirm the topolog- 
ical order numerically, or measure other interesting quantities, 
given a working Monte Carlo simulation. 



A. Correlation functions 

Correlation functions are an obvious starting point. Of 
course, a topological order state has exponentially decaying 
correlations, so this serves primarily as a negative test: we 
check that the system is not a height model in disguise (see 
Sec. lIII Bt . which would have power-law correlations, and that 
it doesn't have long-range order (which can emerge even in 
equal-weighted entropic ensembles, or because the defining 
constraints are too restrictive). Correlations are also of inter- 
est near a critical point where long-range or quasi-long-range 
order emerges. 

In models with vector spins S;, one was accustomed to eval- 
uating the expectation of Sj ■ Sj, or occasionally its second 
moment. It may not be immediately obvious what to mea- 
sure now. One can, of course, simply tabulate frequencies of 
different combinations, e.g. (for the "height difference") how 
often 7o^ r belongs to each conjugacy class. It is preferable, 
though, to reduce the measurements to a single (meaningful) 
number, and the appropriate generalization of the dot product 
is the trace of the matrices in the right group representation. 

Thus we are led to use a character function x( x )> where 
x is any group element; this is always the same within each 
conjugacy class of the group. I divide the actual character by 
the dimension of the representation, so that x( e ) = 1 f° r an Y 
representation, and |x(a;)| < 1 for any element. Presumably, 
the best choice of representation is the one that has the largest 
positive xi 17 ) f° r spins (for a G S). This corresponds concep- 
tually to using a distance metric, within the group Q, counting 
many multiplications by some element of S are needed to take 
you from element to the other one. 

1. Height difference correlation 

In the old "height models" (sketched in Sec II Al l, a natural 
measure of fluctuations was (\h)0) — h(r)\ 2 ). The natural gen- 
eralization of this for the present models with finite (possibly 
non-abelian) groups is 

C h (v) = (x(7(io^r))). (5.1) 

Of course, the product 7(^o^ r ) is independent of which path 
is taken from to r - provided the path does not wrap around 
the periodic boundary conditions. 

As just noted, choosing x(-) so that x{°~) is as close to 
one as possible, provides that Ch(r) does express how fast 
the group element wanders from the identity under repeated 
compositions; that is the choice likeliest to give a monotonic 
decay with distance. If 7(£o-s. r )) is equally likely to be any 
group element - which one expects large r - then it follows 
that C(r) =0. 

2. Spin correlations 

Similarly, we can compute 

Gij = (xi^^aj 1 )). (5.2) 
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B. Defects 

It is easy to augment the simulation to allow a defect pla- 
quette where the plaquette constraint is violated. The same 
(single-site) update rules will work correctly next to the de- 
fect, but they cannot change its position. To make a defect 
mobile, one can add additional update rules specific to the 
defect, by (say) arbitrarily choosing one bond of the plaque- 
tte and changing it to make the plaquette's loop product be 
e (which, of course, the loop product not be e for the pla- 
quette on the other side of that bond, unless that was also a 
defect plaquette and this is the annihilation event.) The sim- 
ulation would normally be run with a constraint or bound on 
the number of defects. 

The idea is to create a pair of defects, by hand, and then 
evaluate expectations depending on them. The first thing to 
measure is the distribution Pd(R) of defect separations R. In 
the case of topological order, we expect deconfinement, mean- 
ing Pd(R) — > const for R > £, where £ is a (not very large) 
correlation length. One can define an effective (entropic) po- 
tential V(R) by 

P d (R) cx exp(V(R)); (5.3) 

physically, V(R) is the difference in entropy due to plac- 
ing the defects near to each other. In the case of a height 
model, V (R) oc In |R|, and dp(R) decays to zero as a power 
law. B 

In fact, since there are various flavors of defect labeled by 
different group elements b, one really needs to write the effec- 
tive potential as 

E = U(b) + U(b r ) + V MV (R) (5.4) 

where b and b' are the respective defect charges, and c is the 
net charge of the combined defect. Here, U(b) and U(b') are 
"core energies" of these respective defects; these, and usually 
the inter-defect potential, are functions only of the conjugacy 
classes of b, b', and/or c. Implicit in the form ( 15.41 is that the 
exponential confinement length probably depends on all of b, 
b', and c. 

Measuring how the effective potential depends on class is 
more physical, since (i) it decides whether a defect is stable 
against decays into other defects (ii) measurements of defect 
behavior (in simulations or in real systems, were any to be 
discovered) might be used to discover the universality class of 
the topological order, if that were not known. I conjecture that 
the dependence on b, b', and c is also described by a character 
function; it would be interesting to see if that can be explored 
analytically in some model. 

Incidentally, since (with the appopriate boundary condi- 
tions) we can have a single defect in our system cell, that gives 
additional opportunities to evaluate e.g. the core energy U(b) 
without the complication of a second defect. 

Finally, if the single-site updates of Sec. lIV Cl are not feasi- 
ble, defects provide a less elaborate alternative update move, 
in place of the cluster update of Sec. IIV Al Namely, we cre- 
ate a pair of defects and allow them to random-walk until they 
annihilate. (However, if their paths differ by a loop around 



the periodic boundary conditions, they may be unable to an- 
nihilate.) Many Monte Carlo schemes I31II32I1 are based on a 
similar process. 

C. Topological sectors 

The tests of topological order outlined up to here have been 
negative; none of them catpures the positive property of topo- 
logical order, which is the degeneracy of topological sectors. 
This can be measured in a classical simulation, if we use a 
(necessarily nonlocal) update which can change sectors, while 
satisfying the detailed balance condition. Either the cluster 
update of Sec. lIV Al or the defect-pair update just outlined in 
Sec.EHwill suffice. 

From the relative fraction of time spent in different topo- 
logical sectors, we can infer a free energy Fl(j x ,j v ), where 
(lx,1y) are the loop products characterizing the sector, and L 
refers to the system size. This is a finite size effect, since (by 
definition of topological order) the difference between sectors 
vanishes in the thermodynamic limit; Fl is expected to decay 
exponentially as a function of L. 

In a similar fashion, if we allow transitions between states 
with and without a defect as part of the dynamics, we can 
evaluate the defect core energy U(b). Of course, Fl(j x ,j v ) 
is very analogous to U(b), since b is a loop product encircling 
the puncture where a defect sits, just as 7^ is from the loop 
product encircling the system. B45I1 

VI. TRANSFER MATRIX AND ANALYTIC APPROACHES 

In a quantum mechanical models with topological order, 
the energy differences between different topological sectors 
decays with system size as cxp(const L), and the correlations 
of a defect pair decay as exp(— Up to now, I have as- 
sumed without justification that this would carry over to the 
present classical models. 

This section finally examines the basis of exponential be- 
havior. I turn here to an analytic treatment using the (practi- 
cally) one-dimensional framework of transfer matrices. First 
of all, this sheds some light on why the finite-size depen- 
dences, as well as the defect-defect interaction, are exponen- 
tially decaying with distance. More specifically, they clar- 
ify the pattern of how sector-weight splittings or defect-pair 
distributions relate to the group's representations and symme- 
tries. 



A. One-dimensional model 

Imagine the most trivial system which can have topologi- 
cal sectors: the one-dimensional version of the discrete-group 
height models. There can be no plaquette constraint. Our en- 
semble simply consists of chains of length L - with periodic 
boundary conditions - having a group element <7j placed on 
each link, the only constraint being that a G S. All {ns) L 
sequences are equally likely. 
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If we let 



j(x) 



a x * a x _i * ... * <j\. 



(6.1) 



then the topological sectors are labeled by j(L). Define the 
(rig x rig dimensional) transfer matrix T in the standard fash- 
ion: let T 7 / >7 be the number of ways to get 7(2; + 1) = 7' 
from 7(2) = 7. Then (T L ) 1:C is the partition function (the 
total number of states) for the sector with j(L) = 7. Note 
that T commutes with permutations that implement the sym- 
metry operations (automorphisms) of the group Q; hence, the 
eigenvalues/eigenvectors of T are classified by the represen- 
tations of the automorphism group (mentioned in Sec. lIII At . 
The transfer matrix has eigenvalues {Aj.} and correspond- 
ing eigenvectors {i>fc, m }; the index m labels each family of 
symmetry-related eigenvectors belonging to the same (degen- 
erate) eigenvalue A&. 

The restricted partition function for topological sector 7 is 



Vk,: 



(6.2) 



Hence, in any sector the overall (entropic) free energy per 
unit length is InAo, where Ao is the largest eigenvalue, 
and L-dependent corrections depend on some larger eigen- 
value A s . For this trivial one-dimensional model, vq — 
(1,1, ...,1,1)/ y/ng. More generally vo must be totally sym- 
metric under all automorphisms of Q, i.e. it belongs to the 
trivial representation. Indeed, v s must also belong to the triv- 
ial representation, since [ufc, m ] e in ( 16.21) is independent of m, 
but 2 m ['ffc, m ] 7 = for any other representation. We let A s 
be next largest (necessarily nondegenerate) eigenvalue of the 
fully symmetric representation, after Ao. 
Hence, 



P( 7 ) l + c 7 (A s /A ) J 



where 



P(e) l + c e (Ai/A ) J 



(6.3) 



(6.4) 



where [wo]g[i>n] e = l/ n S> for this one-dimensional model. It 
follows from (16. 3b that 



In 



P{l) 
P{e) 



(6.5) 



where exp(— l/£i) = |A s /Ao|. Often Ai < 0; in this case, 
we must add a factor (— 1) L on the right-hand-side of ( 16.5b . 
Furthermore, at short L, we may see subdominant terms with 
shorter decay lengths £2 etc., deriving from other eigenvectors 
ofT. 



1. Example 

A useful example is any group Q when S — Q \ e, 
1. every element but the identity is allowed. In this case 



T = (1, 1, 1) <g> (1, 1, 1) - I. Thus A = ng - 1 and 
Ai = A 2 = ... = -1. Thus exp(-l/£i) = l/(n g - 1) and 
the deviations have alternating signs, i.e. the (— 1) L factor is 
needed in ( 16.5b . For the model Z2 x ^(2), the matrix is 



T = 



(l 2 2 2\ 
2 12 2 
2 2 12 

\2 2 2 1/ 



(6.6) 



2. Symmetry-class reduced matrix 

We can classify eigenvectors as "symmetric" or "asymmet- 
ric" according to what representation of the automorphism 
group they transform under. "Symmetric" eigenvectors are 
invariant under group symmetries, while "asymmetric" eigen- 
vectors represent a bias of the probability distribution favoring 
certain local patterns over other (symmetry-related) ones. 

Since the sector probability ratio (16.3b is the same for all 
symmetry-related 7, 1 believe not only v — but also vi must 
be totally symmetric. That affords a considerable simplifica- 
tion, for we can replace T by its projection T onto the group 
element symmetry classes. (Such a class consists of elements 
that map to each under under some automorphism, so these 
are at least as large as the conjugacy classes.) Whereas the 
dimension of T was the number of group elements ng, the 
dimension of T is the number of group classes: Tji tells the 
number of times that 177 belongs to class j, if 7 belongs to 
class i and a runs over all rig elements in 5. 

For example, in the case of the group A5, the matrix is re- 
duced from 60 x 60 to 4 x 4, with entries for elements of 
order one (identity), two, three, and five. (There are two con- 
jugacy classes with order five, but they are equivalent by an 
outer automorphism.) For the model A$(3), we get 



f 1 o\ 

4 6 5 

20 8 7 5 

\ 8 6 10/ 



(6.7) 



This matrix is similar to a symmetric matrix L) _1 / 2 T£) 1 / 2 ; 
here D = diag(l, 15, 20, 24) for this group, or in general is 
the diagonal matrix with entries being the count of each class. 



B. Sector probabilities in two dimensions? 

A two-dimensional finite-group height model is also de- 
scribed by a transfer matrix T. However, now the vector that 
T acts on represents all possible path products j x , v taken to a 
point (x,y), and thus is (ng) w dimensional, where W is the 
width of the strip (in the y direction; iteration still runs in the x 
direction). We must replace c 7 — > c^(W) and £1 — > £(W) in 
Eq. ( |6.5b . Conceivably £(W) — > as W — > 00, as is very well 
known in gapless systems, so the form of exp(— L/£(W)) 
does not prove exponential decay in d = 2. 
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Nevertheless, we can make a plausible guess to obtain a 
fitting form for comparison with numerics. Since all correla- 
tions are expected to be rapidly decaying, a strip of width W 
is like W/wq independent, one-dimensional strips of width 
Wq in parallel. But all these strips are constrained to have the 
same, or equivalent, sector label 7.) The consequence is that 



P(e) 



l + c 7 (A 1 /A ) i 

l + C e (Al/A ) 2 



W/wo 



In 



P(e) 



(C 7 - C e )We 



-L/t 



(6.8) 



(6.9) 



c-y/wo and £ w £1 independent 



in place of ( 16.51 ). with C 7 
of W. 

My chief motivation for introducing the transfer-matrix for- 
malism is separate from such guesses about the W scaling, 
and is much better founded. Namely, the eigenvectors for 
the corrections to P(j) are representations of the automor- 
phism group. Furthermore, which representation goes with 
the longest correlations is probably the same as in the one- 
dimensional case. What really matters here is that our choice 
of a selected set S defines a sort of metric on Q: the distance 
from g to g' is the number of times you need to multiply by 
an element of S to get from g to g'. Then, the first nontriv- 
ial eigenvector v\ is the mode that is slowest varying on Q 
according to this metric (apart from vq which is uniform). 

The one-dimensional correlation length £ 1 can be computed 
for any combination of Q and S and can serve as another "fig- 
ure of merit" for a group. That is, in light of the previous para- 
graph's argument, it should be roughly related to the true sec- 
tor probability decay length £ for the two-dimensional model 
(and likely related to the defect-defect decay length as well). 



C. Defect separations and W — 2 transfer matrix 

Whereas the one-dimensional model already seems to cap- 
ture the essence of how sector probabilities depend on system 
size and sector label, it does not admit topological defects and 
hence sheds no light on the parallel question of how p(R) for 
a defect pair decays with separation or depends on the respec- 
tive defect charges. 

Clearly, p(R) must be associated somehow with the eigen- 
vectors and eigenvalues of the two-dimensional transfer ma- 
trix, since all possible correlation information is expressed in 
it. But it is not self-evident just what kind of distortion of the 
ensemble is being propagated, or what sort of subdominant 
eigenvector: the symmetric kind (which governed the sector 
probabilities) or the asymmetric kind. 

I will work out here a toy calculation, again using a transfer 
matrix, of the correlation decay due to asymmetric eigenvec- 
tors. I believe they are the ones that matter for the case of 
an abelian group. In that case, the charge of a defect is a 
particular element: the loop product around the defect gives 
that same result, no matter how big the loop, and only another 



defect with the inverse of that charge can cancel it. In the 
non-abelian case, however, these properties would seem to be 
defined only modulo conjugacy classes. Therefore, the picture 
presented here is only asserted to go with abelian groups. 

The simplest property that could influence or be influenced 
by a defect's presence is the correlation of two adjacent spins 
on the same plaquette, i.e. sitting on bonds that make a 90° 
angle. A simple example is the model Z2 x Z2(2)sq, in which 
ns = 3 elements are allowed - all except the identity. These 
elements are {a,b,c}. Consider a plaquette with the spins 
on two edges specified and the remaining two spins to be as- 
signed (there are 3 2 unconstrained ways to do so). When ad- 
jacent edges on a plaquette have the same element, there are 
three ways to satisfy the plaquette constraint, but only two 
ways if the given adjacent edges are different. On the other 
hand, if we want to make a defect plaquette, there are six ways 
when the given spins are the same but seven ways when they 
are different. 

Let's set up a W = 2 strip, the narrowest kind that can 
capture defect correlations. This transfer matrix, unlike the 
previous one, refers to the actual spin configurations in each 
vertical pair of bonds; we add up all the possible horizontal 
bonds. I assume the upper row of plaquettes are constrained 
to be have identity product around the plaquette. Plaquettes in 
the lower row can have any product - defects are permitted - 
with a weight 9q for the identity or 9 a ,9b, C for the respective 
defect charges a, b, c. We imagine the limit in which 9 a ^ c are 
small and ask for the corresponding defect correlations. 

Although T has 3 2 x 3 2 = 81 elements, in fact there are 
only ten distinct kinds by symmetry, as given in Table |IV| 
"no." represents the number of times each kind occurs in the 
matrix. The factors 9q w 1 and 9 a <C 1 are omitted in the ta- 
ble. To compute the matrix elements, note that when o\ = <r[ 
in the upper plaquette, the central horizontal bond may be any 
element [three possibilities] but if o\ ^ a[, the central bond 
may be only a\ or <t[ [two possibilities]. There are always 
three possibilities for the lower horizontal bond, so the table's 
rows add up to 9 or 6 depending whether or not o\ = <j[ . 

The probability to find a defect of charge f3' at separation 
R, given there is a defect of charge (3 at the origin, is then 



p(R) 



Tr 



L-R-l T (/3') [ T (0)]fl-l T (/3)^ 



Tr([T(°)] L - 1 T(/ 3 ) 



(6.10) 



For a large power M, we can replace [T^°>] M — > (A ) M i>o <8> 
i>o + (Ai) M t>i (£> V\. Here v and v\ are the eigenvectors be- 
longing to the maximum and next-largest eigenvalues of T'"'. 
Assuming L ^> i? ^> 1, we get 

Aj- fl - 1 i: t „, 1 EJ''o,T^, m )^, m r^o)A«- m 1 



p(R) = 



A^KT^o) 



where 



MP) 



A 



(6.11) 
(6.12) 

(6.13) 
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TABLE IV: Example for group Zi x Zi{T): Transfer matrix elements 

ji(fS) 
°'i' fJ o ;0 ' 1,cro 



no. 




<T ) (a[,a' ) 


j,(0) 


Jl((t) rp(b) rp(c) 


3 


(aa) 


(aa) 


3 


2 


2 2 


12 


(aa) 


(ab) 


2 


2 


2 3 


12 


(aa) 


(ba) 


2 


1 


1 2 


6 


(aa) 


(bb) 


2 


1 


1 2 


12 


(aa) 


(be) 


1 


2 


2 1 


6 


(ab) 


(ab) 


3 


2 


2 2 


6 


(ab) 


(ba) 


2 


1 


1 2 


6 


(ab) 


(ac) 


2 


3 


2 2 


6 


(ba) 


(ca) 


2 


2 


1 1 


12 


(ab) 


(ca) 


1 


2 


1 2 



- remember (uq, T^vq) = Aq] — and 



/A y/2 Kr ( ^o) 



(6.14) 



Please remember, the eigenvector called «i here is asym- 
metric, and is thus not the same as the symmetric eigenvector 
called v s in Sec. lVI Al We see that asymptotically, 



\np(R) oc c(/3)c(/3')e 



-R/i 



(6.15) 



Notice first that the decay length £ is independent of the defect 
charges, but different defect charges have different projections 
onto this eigenmode. As is clear from the derivation, a more 
general form could be written, including subdominant contri- 
butions: 



D. Other approaches to p(R) in d — 2 

I conjecture there is an alternative approach which is more 
congenial to d = 2. Namely, in the vicinity of a defect, the 
probabilities of local patterns have small deviations from the 
bulk values, which could be represented by operators Ok (r) 
and small conjugate fields /ifc(r). That is, adding a Hamil- 
tonian ^ r ft,fe(r)Ofc(r) (in the absence of the defect) would 
perturb the ensemble the same way as the defect does. Note 
that the operator "Ok (r)" is schematic, in the sense that such 
operators probably involve two spins at different r (in light of 
the same logic laid out in the first paragraphs of this subsec- 
tion). Then possibly some sort of mean-field approximation 
produces a difference equation for /ifc(r), the discrete ana- 
log of Poisson's equation V 2 /i/c(r) = /ifc(r)/^, which has 
solutions ~ e~ R ^ k /R. In this approach, we have a sort of 
small parameter in that the influence of a perturbation decays 
as e~ R ^ k , which becomes arbitrary small at sufficiently large 
R. We can therefore rely on linear response in that regime. 

One can conceive additional approaches to p(R) which 
depend on a genuine small parameter; the difficulty is that 
the actual model families defined in this paper are far from 
that limit. For example, one could expand around the pure 
gauge theory: in place of the spin constraint, there would 
be no constraint but configurations would have a statistical 
weight exp(A^ r r < u ( cr ( r i r '))' where u(a) would penalize 
all a ^ S. In the limit A — > 0, we have a pure gauge the- 
ory in which all correlation lengths are zero, so hopefully 
£fc would scale as a power of A. The models under consid- 
eration are, unfortunately, the case A = oo. Still, since the 
topological phases are like the pure gauge models at large 
scales, they should be adiabatically connected and hence this 
approach should be qualitatively valid. 



Inp(R) cx c k (P)c k (P')e- R '^ (6.16) 

k 

where £i > £2 > The later terms could be important 
corrections to include in fits at short R, particularly when the 
smaller ^'s happen to be associated with larger coefficients 
Ck{/3). Also, if Ci(/3) = for certain defects, their asymptotic 
interaction gets carried by the first mode that has nonzero pro- 
jections onto both defects. 

The formulas basically apply to any width of strip. (If de- 
fects are allowed in more than one horizontal row of plaque- 
ttes, then the defect distribution is no longer a function just 
of R but also of the two y coordinates; the only modification 
necessary is that -» T^' v \ labeled not only by the de- 
fect's flavor but by its y coordinate.) I would speculate that 
the W = 3 strip, with defects in the central row, may be a 
good approximation in practice, although of course there is 
no control parameter to make small. The basis for this is sim- 
ply the notion that, when we have rapid exponential decays, 
these are associated in the ensemble with strings connecting 
the defects; any influence carried by a less direct chain would 
be exponentially smaller in correspondence with the longer 
length. 



VII. DISCUSSION 

I have put forward the notion of purely classical topological 
order, defined by an ergodicity breaking into sectors depen- 
dent on the topology, and not distinguishable by thermody- 
namic expectations of any local operator. A family of explicit 
models has been described, along with a suitable Monte Carlo 
technique, and criteria were suggested to pick out the most 
promising cases (in having a nontrivial and updatable ensem- 
ble of allowed states). 

A framework was set up (Sec. Ill Al ) to define models with 
three variable attributes: which group, which class(es) out of 
the group to selected as "spin" variables, and which lattice 
to place the model on. These are characterized by parame- 
ters - the sizes ng and ns of the group and the spin subset, 
the coordination number z of the lattice and Zd of its dual 
- which entered crude formulas that estimate the entropy of 
the model and its updatability under single-site Monte Carlo 
moves (Sec. IIIICl and IIV Bl ). which are the only actual cal- 
culations in the paper. Groups with normal subgroups tend to 
be "less nonabelian ", thus perhaps less attractive (Sec. lIIIBl . 
Although topological order superficially would appear to be 
intrinsically featureless, there is sufficient richness of measur- 
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able functions when one considers the dependence of free en- 
ergy on topological indices - finite size dependence on sector 
or finite distance dependence on defect separations (Sec. Ml. 

In trying to connect the classical picture to the quantum 
theory of topological order, it is intriguing that a given two 
defect charges (see Sec. Ill Q can combine in more than one 
way, in a classical non-abelian model, reminiscent of fusion 
rules in a quantum model. If further investigation finds that the 
sector counting gives the same degeneracies in the classical as 
in the quantum case, one would conclude that this is one of 
the shared properties, not an intrinsically quantum one. 

The physical manifestations of classical topological order 
and/or of non-abelianness are less striking, perhaps, than for 
the quantum case. Most prominent is the behavior of topolog- 
ical defects. Topological order implies deconfinement in the 
classical model for nonabelian and abelian cases alike. Non- 
abelianness (of the group) changes the rules for addition of 
defect charges, and braiding has physical consequences, even 
though there are no Berry phases in a classical model. (It 
must be noted, however, that the same behaviors are seen in 
non-abelian defects of ordinary long-range order [28] - they 
are not inherent to topological order.) 

Degeneracies of different topological sectors, the defining 
property of topological order, work differently in the non- 
abelian than in the abelian case: for example, there are far 
distinct fewer sectors in the non-abelian case (Sec. Ill Bl i. 



A. Quantum mechanics 

Several central concepts of topological order are inherently 
quantum-mechanical and thus have no counterpart in classi- 
cal topological order. They are mainly related to phases in 
wavefunctions and braiding of worldlines in 2+1 dimensions, 
namely anyon and mutual statistics. Most real or imagined ex- 
periments relating to topological excitations have involved in- 
terference phenomena (e.g. tunneling in various geometries of 
quantum Hall fluids) and thus probe the quantum-mechanical 
aspects of topological order. 

Another feature missing in the classical models is the dual 
defect or quasiparticle (such as the vison IB9IP . which is a 
distortion of the phase factors in the many-body wavefunction. 

A final attribute of topological orders is the nontrivial 
counting statistics of the excited states made by several quasi- 
particles, which is quantum mechanical in that it concerns the 
linear dimension of a Hilbert space. One cannot rule out the 
appearance of similar concepts in classical stat mech - there, 
too, the partition function contains combinatorial factors for 
the placement of defects, after the other degrees of freedome 
have been integrated out. However, I am not aware of a clas- 
sical situation in which such a nontrivial counting actually 
emerges. 

B. Construction of quantum models? 

Any of the classical height models with topological order 
may be converted into a simular quantum model if we can 



endow it "flipping" move, just as classical dimer (and other) 
models get converted into quantum dimer models using the 
Rokhsar-Kivelson (RK) prescription l37tl . A barrier to this 
is that the only generally guaranteed "flip" move is a cluster 
update, as explained in Sec. lIV Al 

Fortunately, whenever the single-site update (Sec. IIVCI ) 
suffices, we can define a quantum model with a simple "flip- 
ping" term in the Hamiltonian, usually parametrized by an 
amplitude t, as well as a "potential" term of strength V = t 
that penalizes each flippable place. At the RK point V = t, 
the ground state wavefunction is a superposition of all config- 
urations in the same topological sector, with the same (equal) 
weighting as in the classical ensemble, and (mutually inac- 
cessible) topological sectors are trivially degenerate. One is 
also free to set V — - obtaining a simpler model in which 
flippable sites are so favored that an ordered state is likely to 
be the outcome - or to vary V/t with the hope of crossing a 
phase transition. 

The above recipe is incomplete, in that there are many pos- 
sible choices of update (labeled by the multiplier r of Sec.lIVI), 
and presumably all or many should be included in "flipping" 
term of the quantum Hamiltonian, which requires a prescrip- 
tion for the relative magnitudes of coefficient to put for each 
class of t, as well as the relative phase factors. Presumably, a 
proper choice is taking the same phase factors for every term, 
i.e. the Hamiltonian transforms by the fully symmetric (triv- 
ial) representation of the automorphism group of Q. Alterna- 
tively, in lucky cases, one might select a site-dependent pat- 
tern of Tj's so as to link the group symmetry to the lattice sym- 
metry, in the spirit of Kitaev's honeycomb model |4|]. Another 
option is to include a second, quantum fluctuating field of t's 
which are used for the update. If the t's are derived from a 
second set of "spins" also having the gauge-like structure of 
a finite-group height model, we might be able to realize dual 
("magnetic") defects having mutual statistics with the er-spin 
type ("electric") defects described in this paper. 



C. Possible simulations 

As laid out in Sec. |Vj several quantities e.g. correlation 
functions can be measured in classical non-abelian height 
models as a test-bed, whereas the analogous calculation might 
be very challenging computationally in a quantum mechanical 
model. Of course, the answers need not be the same, but the 
questions may be much clearer once the classical results are 
in hand. First, one can create defect pairs and evaluate the his- 
togram of their separations, which will reveal whether or not 
they are deconfined. Secondly, one can evaluate the probabil- 
ities of different topological sectors, which is the direct test of 
topological order. 

Furthermore, if we generalize the models to include classi- 
cal Hamiltonians (so as to weight configurations according to 
the Boltzmann distribution), phase transitions can be studied. 
Just as a standard height model may have a "smooth" phase, 
in which one or more height components becomes locked, it 
seems conceivable that a discrete-group height model might 
have a phase in which the loop products can take values in a 
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subgroup Q' C Q . If so, one might encounter critical points 
separating different topological phases, and characterize the 
critical exponents. 

D. Dilution and effective interactions of local degrees of 
freedom? 

A classical model might be a helpful too for investigating 
the consequences of dilution disorder in a model with topolog- 
ical order. Each diluted site or plaquette is like a very small 
hole cut in the system, thereby increasing the genus and the 
number of topological sectors. If the hole were big, these sec- 
tors would be truly degenerate (by the definition of topological 
order), however this degeneracy is broken since the holes are 
small. 

The values of 7* on each dilution site are local pseudospins, 
which are expected to have (exponentially decaying) interac- 
tions mediated by the fluid in between them, like the emergent 
spin- 1/2 degrees of freedom in diluted spin-1 antiferromag- 
netic chains |f33tl . The ground state of such a system could be 
constructed by a renormalization group that iteratively com- 



bines the most strongly coupled pair of pseudospins into a sin- 
gle effective pseudospin, as was originally done for the (expo- 
nentially decaying) antiferromagnetic coupling of the charge- 
bound electron spins in P-doped Si [36]. 

In the non-abelian case, at least, we do not know for sure 
whether these interactions lead to an inert singlet phase (as in 
the antiferromagnet) or could give a state with some kind of 
order among the pseudospins. Thus the system with defects 
would not be a topological liquid, and this would be a novel 
scenario of how order can emerge due to disorder. l38ll 
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